# Summary: Creates Appendix I

######################################################
#-------------------Appendix I ----------------------#
######################################################

options(scipen=999)
rm(list = ls())

# Load libraries
library(rddensity)
library(rdrobust)
library(rdpower)
library(sandwich)
library(stargazer)
library(readxl)
library(lmtest)

# Set WD
# setwd('~/Dataverse/')

# Source functions
source('./Code/fun_analysis.R')

# Load Data
load(file = "./Data/d.RData")
load(file = "./Data/d_balance.RData")

# Combine datasets
d <- merge(d[, c('codeyear', 'NUMERO_PARTIDO', 'position',
		'SIGLA_PARTIDO', 'run_var', 'eleito', 'CODIGO_SEXO', 'diff_cand', 'mayorFem')]
	, d_balance, by = c('codeyear', 'NUMERO_PARTIDO'), all.x = TRUE)


# Variables removing laranjas
d$avg_DF_nolar <- ((
	(d$total_votes_party_DF-d$total_votes_party_laranja_DF)/d$total_votes_district_DF) * 100) /
	(d$n_candidates_DF - d$n_candidates_laranja_DF)

d$avg_fem_DF_nolar <- ((
	(d$total_votes_party_fem_DF-d$total_votes_party_fem_laranja_DF)/d$total_votes_district_DF) * 100) /
	(d$n_candidates_fem_DF - d$n_candidates_fem_laranja_DF)

d$avg_male_DF_nolar <- ((
	(d$total_votes_party_DF-(d$total_votes_party_fem_DF -d$total_votes_party_fem_laranja_DF)- d$total_votes_party_laranja_DF)/d$total_votes_district_DF) * 100)/
	(d$n_candidates_DF - (d$n_candidates_fem_DF - d$n_candidates_fem_laranja_DF) - d$n_candidates_laranja_DF)

# Women
d$n_cand_w <- d$n_candidates_fem_DF - d$n_candidates_fem_laranja_DF
d$n_cand_m <- (d$n_candidates_DF - (d$n_candidates_fem_DF - d$n_candidates_fem_laranja_DF) - d$n_candidates_laranja_DF)

# load(file = "Data/d.RData")
d$nonZeroDF <- NULL
zerosDF <- aggregate(d$n_cand_w>0 & d$n_cand_m>0, list(d$codeyear), sum)
names(zerosDF) <- c('codeyear', 'nonZeroDF')

# Add to dataset
d <- merge(d, zerosDF, by = 'codeyear')

# Keep only observations from municipalities where both parties had male and female candidates
d1 <- subset(d, nonZeroDF == 2)

# Get codeyear for where a centralized and decentralized party won
centParties <- d1$codeyear[(d1$SIGLA_PARTIDO %in% c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT")) & d1$eleito == 1]
decentParties <- d1$codeyear[!(d1$SIGLA_PARTIDO %in% c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT")) & d1$eleito == 1]

#################################################################
####### Only Centralized Parties are able to deliver votes ######
#################################################################

# FIGURE I.60
pdf('./AppendixIResults/figI60.pdf')
# All Parties
genPlot(mydata = subset(d1), 
	DV = 'avg_DF_nolar', 
	dvName = 'Federal Deputy', 
	vs = T, 
	clusterVar = subset(d1)$codeyear)
# Centralized
genPlot(mydata = subset(d1, codeyear %in% centParties), 
	DV = 'avg_DF_nolar', 
	dvName = 'Federal Deputy', 
	vs = T, 
	clusterVar = subset(d1, codeyear %in% centParties)$codeyear)
# Descentralized
genPlot(mydata = subset(d1, codeyear %in% decentParties), 
	DV = 'avg_DF_nolar', 
	dvName = 'Federal Deputy', 
	vs = T, 
	clusterVar = subset(d1, codeyear %in% decentParties)$codeyear)
dev.off()

# Main Models: h = optimal 
out.all.opt <- rdrobust(y = subset(d1)$avg_DF_nolar
		, x = subset(d1)$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1)$codeyear
		, all = TRUE)
out.cent.opt <- rdrobust(y = subset(d1, codeyear %in% centParties)$avg_DF_nolar
		, x = subset(d1, codeyear %in% centParties)$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties)$codeyear
		, all = TRUE)
out.decent.opt <- rdrobust(y = subset(d1, codeyear %in% decentParties)$avg_DF_nolar
		, x = subset(d1, codeyear %in% decentParties)$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties)$codeyear
		, all = TRUE)
# Test Difference:
zOut <- zDif(beta1 = out.cent.opt$coef[3,], beta2 = out.decent.opt$coef[3,]
	, se1 = out.cent.opt$se[3,], se2 = out.decent.opt$se[3,])
zOut
(1 - pnorm(zOut$zscore)) * 2

#### Results by Period
d1$period <- 'weak'
d1$period[d1$ANO_ELEICAO > 2007] <- 'strong'
d1$period[d1$ANO_ELEICAO > 2015] <- 'moderate'

# Main Models: h = optimal 
out.cent.opt.weak <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'weak')$avg_DF_nolar
		, x = subset(d1, codeyear %in% centParties & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'weak')$codeyear
		, all = TRUE)
out.cent.opt.strong <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'strong')$avg_DF_nolar
		, x = subset(d1, codeyear %in% centParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'strong')$codeyear
		, all = TRUE)
out.cent.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'moderate')$avg_DF_nolar
		, x = subset(d1, codeyear %in% centParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'moderate')$codeyear
		, all = TRUE)
out.decent.opt.weak <- rdrobust(y = subset(d1, codeyear %in% decentParties & period == 'weak')$avg_DF_nolar
		, x = subset(d1, codeyear %in% decentParties & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties & period == 'weak')$codeyear
		, all = TRUE)
out.decent.opt.strong <- rdrobust(y = subset(d1, codeyear %in% decentParties & period == 'strong')$avg_DF_nolar
		, x = subset(d1, codeyear %in% decentParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties & period == 'strong')$codeyear
		, all = TRUE)
out.decent.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% decentParties & period == 'moderate')$avg_DF_nolar
		, x = subset(d1, codeyear %in% decentParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties & period == 'moderate')$codeyear
		, all = TRUE)

# Get point estimates, se, and ci for men
betas_cen <- c(out.cent.opt.weak$coef[3,], out.cent.opt.strong$coef[3,]
	, out.cent.opt.moderate$coef[3,])
se_cen <- c(out.cent.opt.weak$se[3,], out.cent.opt.strong$se[3,]
	, out.cent.opt.moderate$se[3,])
lb_cen <- c(out.cent.opt.weak$ci[3,1], out.cent.opt.strong$ci[3,1]
	, out.cent.opt.moderate$ci[3,1])
ub_cen <- c(out.cent.opt.weak$ci[3,2], out.cent.opt.strong$ci[3,2]
	, out.cent.opt.moderate$ci[3,2])

# Get point estimates, se, and ci for women
betas_decen <- c(out.decent.opt.weak$coef[3,], out.decent.opt.strong$coef[3,]
	, out.decent.opt.moderate$coef[3,])
se_decen <- c(out.decent.opt.weak$se[3,], out.decent.opt.strong$se[3,]
	, out.decent.opt.moderate$se[3,])
lb_decen <- c(out.decent.opt.weak$ci[3,1], out.decent.opt.strong$ci[3,1]
	, out.decent.opt.moderate$ci[3,1])
ub_decen <- c(out.decent.opt.weak$ci[3,2], out.decent.opt.strong$ci[3,2]
	, out.decent.opt.moderate$ci[3,2])

# Calculate Differences 
zOutDecent <- zDif(beta1 = out.decent.opt.strong$coef[3,], beta2 = out.decent.opt.weak$coef[3,]
	, se1 = out.decent.opt.strong$se[3,], se2 = out.decent.opt.weak$se[3,])
zOutCent <- zDif(beta1 = out.cent.opt.strong$coef[3,], beta2 = out.cent.opt.weak$coef[3,]
	, se1 = out.cent.opt.strong$se[3,], se2 = out.cent.opt.weak$se[3,])


# Start the plot:
# FIGURE I.61
pdf('./AppendixIResults/figI61.pdf', width = 10)
par(mar = c(5, 5, 2, 2))
plot(x = NULL, y = NULL, axes = FALSE, 
	cex = 2, xlim = c(0.75, 2.25), ylim = c(-1, 1),
	ylab = "RD Effect of Co-Partisan Mayor", col = c('black', 'grey46'), 
	pch = c(20, 17),
	xlab = "Party Type", 
	cex.lab = 2)
sapply(1:2, function(i) lines(x = c(c(0.9, 1.1)[i], c(0.9, 1.1)[i]), 
	y = c(lb_decen[i], ub_decen[i]), col = c('black', 'gray46')[i]))
sapply(1:2, function(i) lines(x = c(c(1.9, 2.1)[i], c(1.9, 2.1)[i]), 
	y = c(lb_cen[i], ub_cen[i]), col = c('black', 'gray46')[i]))
axis(1, at = seq(1, 2, 1), label = c('Decentralized', 'Centralized'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_decen[1], y1 = betas_decen[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_decen[1], y1 = betas_decen[2], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_decen[2], y1 = betas_decen[2], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_cen[1], y1 = betas_cen[1], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_cen[1], y1 = betas_cen[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_cen[2], y1 = betas_cen[2], lty = 2)

# Add points for centralized
points(x = c(1.9, 2.1), y = betas_cen[1:2], pch = c(20, 17), col = c('black', 'grey46'), cex = 2)
# Add points for decentralized
points(x = c(0.9, 1.1), y = betas_decen[1:2], pch = c(20, 17), col = c('black', 'grey46'), cex = 2)


# Add information about difference
text(round(zOutDecent$difference, 3), y = betas_decen[1] + 0.7, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutDecent$zscore > 0, (1 - pnorm(zOutDecent$zscore)) * 2, pnorm(zOutDecent$zscore) * 2), 3), ')')
, y = betas_decen[1] + 0.5, x = 1, cex = 1.5)

text(round(zOutCent$difference, 3), y = betas_cen[1] + 0.75, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutCent$zscore > 0, (1 - pnorm(zOutCent$zscore)) * 2, pnorm(zOutCent$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutCent$zscore > 0, (1 - pnorm(zOutCent$zscore)) * 2, pnorm(zOutCent$zscore) * 2), 3)), ')')
, y = betas_cen[1] + 0.6, x = 2, cex = 1.5)

# Add legend
legend('bottomleft', 
	title = 'Electoral Rule Regime:',
	legend = c('Strong (2010-2014)', 'Weak (1998-2006)'), 
	col = c('gray46', 'black'), 
	pch = c(17, 20), 
	lty = c(1, 1),
	bty = 'n',
	cex = 1.5)
dev.off()


################################################################################
####### H1: male candidates > female candidates By Period, only centralized ######
################################################################################

d1$period <- 'weak'
d1$period[d1$ANO_ELEICAO > 2007] <- 'strong'
d1$period[d1$ANO_ELEICAO > 2015] <- 'moderate'

# Weak
out.cent.men.opt.weak <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'weak')$avg_male_DF_nolar
		, x = subset(d1, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'weak')$codeyear
		, all = TRUE)
out.cent.women.opt.weak <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'weak')$avg_fem_DF_nolar
		, x = subset(d1, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'weak')$codeyear
		, all = TRUE)

# Strong
out.cent.men.opt.strong <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'strong')$avg_male_DF_nolar
		, x = subset(d1, codeyear %in% centParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'strong')$codeyear
		, all = TRUE)
out.cent.women.opt.strong <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'strong')$avg_fem_DF_nolar
		, x = subset(d1, codeyear %in% centParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'strong')$codeyear
		, all = TRUE)

# Moderate
out.cent.men.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'moderate')$avg_male_DF_nolar
		, x = subset(d1, codeyear %in% centParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'moderate')$codeyear
		, all = TRUE)
out.cent.women.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'moderate')$avg_fem_DF_nolar
		, x = subset(d1, codeyear %in% centParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'moderate')$codeyear
		, all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.weak$coef[3,], out.cent.men.opt.strong$coef[3,]
	, out.cent.men.opt.moderate$coef[3,])
se_men <- c(out.cent.men.opt.weak$se[3,], out.cent.men.opt.strong$se[3,]
	, out.cent.men.opt.moderate$se[3,])
lb_men <- c(out.cent.men.opt.weak$ci[3,1], out.cent.men.opt.strong$ci[3,1]
	, out.cent.men.opt.moderate$ci[3,1])
ub_men <- c(out.cent.men.opt.weak$ci[3,2], out.cent.men.opt.strong$ci[3,2]
	, out.cent.men.opt.moderate$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.weak$coef[3,], out.cent.women.opt.strong$coef[3,]
	, out.cent.women.opt.moderate$coef[3,])
se_women <- c(out.cent.women.opt.weak$se[3,], out.cent.women.opt.strong$se[3,]
	, out.cent.women.opt.moderate$se[3,])
lb_women <- c(out.cent.women.opt.weak$ci[3,1], out.cent.women.opt.strong$ci[3,1]
	, out.cent.women.opt.moderate$ci[3,1])
ub_women <- c(out.cent.women.opt.weak$ci[3,2], out.cent.women.opt.strong$ci[3,2]
	, out.cent.women.opt.moderate$ci[3,2])

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.cent.men.opt.weak$coef[3,], beta2 = out.cent.women.opt.weak$coef[3,]
	, se1 = out.cent.men.opt.weak$se[3,], se2 = out.cent.women.opt.weak$se[3,])
zOutStrong <- zDif(beta1 = out.cent.men.opt.strong$coef[3,], beta2 = out.cent.women.opt.strong$coef[3,]
	, se1 = out.cent.men.opt.strong$se[3,], se2 = out.cent.women.opt.strong$se[3,])
zOutModerate <- zDif(beta1 = out.cent.men.opt.moderate$coef[3,], beta2 = out.cent.women.opt.moderate$coef[3,]
	, se1 = out.cent.men.opt.moderate$se[3,], se2 = out.cent.women.opt.moderate$se[3,])

# Start the plot:
# FIGURE I.62
pdf('./AppendixIResults/figI62.pdf', width = 14)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
	pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-3, 3),
	ylab = "RD Effect of Co-Partisan Mayor",
	xlab = "Electoral Rule Regime", 
	cex.lab = 2)
sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
	y = c(lb_men[i], ub_men[i])))
sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
	y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# Moderate
segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

# Add points for women's coefficients
points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOutWeak$difference, 3), y = betas_men[1] + 0.7, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
, y = betas_men[1] + 0.5, x = 1, cex = 1.5)

text(round(zOutStrong$difference, 3), y = betas_men[2] + 0.75, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
, y = betas_men[2] + 0.45, x = 2, cex = 1.5)

text(round(zOutModerate$difference, 3), y = betas_men[3] + 0.5, x = 3, cex = 1.5)
text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
, y = betas_men[3] + 0.2, x = 3, cex = 1.5)

# Add legend
legend(x = 0.7, y = -1.5, 
	legend = c('Women', 'Men'), 
	col = c('gray46', 'black'), 
	pch = c(17, 20), 
	lty = c(1, 1),
	bty = 'n',
	cex = 1.5)
dev.off()


#####################################################################################
####### H2: male brokers > female brokers on Male VS, only centralized parties#######
#####################################################################################

# We only want mixed-gender races
d$female <- as.numeric(d$CODIGO_SEXO == 4)
mixed <- aggregate(d$female, list(d$codeyear), sum)
names(mixed) <- c('codeyear', 'mixed')
mixed$mixed <- ifelse(mixed$mixed == 1, 1, 0)
d <- merge(d, mixed, by = 'codeyear')
# Subset to Mixed Races and parties ran both male and female candidates
d2 <- subset(d, mixed == 1 & nonZeroDF == 2)
# Create running variable
d2$run_var_fem <- ifelse(d2$position == 1 & d2$CODIGO_SEXO == 4
	, d2$diff_cand * 1, d2$diff_cand * -1)
d2$run_var_fem[d2$mayorFem == 0] <- NA
d2$run_var_male <- ifelse(d2$position == 1 & d2$CODIGO_SEXO == 2
	, d2$diff_cand * 1, d2$diff_cand * -1)
d2$run_var_male[d2$mayorFem == 1] <- NA
table(is.na(d2$run_var_male))
table(is.na(d2$run_var_fem))
table((d2$run_var_male)>0)
table((d2$run_var_fem)>0)
# Add period
d2$period <- 'weak'
d2$period[d2$ANO_ELEICAO > 2007] <- 'strong'
d2$period[d2$ANO_ELEICAO > 2015] <- 'moderate'


##############################
######### mayor = men ########
##############################
# Weak
out.h2.men.men.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$avg_male_DF_nolar
		, x = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$run_var_male
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$codeyear
		, all = TRUE)

out.h2.men.women.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$avg_fem_DF_nolar
		, x = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$run_var_male
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$codeyear
		, all = TRUE)

# Difference 
outZ <- zDif(beta1 = out.h2.men.men.non.opt$coef[3,], beta2 = out.h2.men.women.non.opt$coef[3,]
	, se1 = out.h2.men.men.non.opt$se[3,], se2 = out.h2.men.women.non.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2
summary(out.h2.men.men.non.opt)
summary(out.h2.men.women.non.opt)

# Strong
out.h2.men.men.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$avg_male_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$run_var_male
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$codeyear
		, all = TRUE)
out.h2.men.women.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$avg_fem_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$run_var_male
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$codeyear
		, all = TRUE)
summary(out.h2.men.men.with.opt)
summary(out.h2.men.women.with.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.men.men.with.opt$coef[3,], beta2 = out.h2.men.women.with.opt$coef[3,]
	, se1 = out.h2.men.men.with.opt$se[3,], se2 = out.h2.men.women.with.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2


# Moderate
out.h2.men.men.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$avg_male_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$run_var_male
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$codeyear
		, all = TRUE)
out.h2.men.women.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$avg_fem_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$run_var_male
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$codeyear
		, all = TRUE)
summary(out.h2.men.men.withLeg.opt)
summary(out.h2.men.women.withLeg.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.men.men.withLeg.opt$coef[3,], beta2 = out.h2.men.women.withLeg.opt$coef[3,]
	, se1 = out.h2.men.men.withLeg.opt$se[3,], se2 = out.h2.men.women.withLeg.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2


# Get point estimates, se, and ci for men
betas_men <- c(out.h2.men.men.non.opt$coef[3,], out.h2.men.men.with.opt$coef[3,]
	, out.h2.men.men.withLeg.opt$coef[3,])
se_men <- c(out.h2.men.men.non.opt$se[3,], out.h2.men.men.with.opt$se[3,]
	, out.h2.men.men.withLeg.opt$se[3,])
lb_men <- c(out.h2.men.men.non.opt$ci[3,1], out.h2.men.men.with.opt$ci[3,1]
	, out.h2.men.men.withLeg.opt$ci[3,1])
ub_men <- c(out.h2.men.men.non.opt$ci[3,2], out.h2.men.men.with.opt$ci[3,2]
	, out.h2.men.men.withLeg.opt$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.h2.men.women.non.opt$coef[3,], out.h2.men.women.with.opt$coef[3,]
	, out.h2.men.women.withLeg.opt$coef[3,])
se_women <- c(out.h2.men.women.non.opt$se[3,], out.h2.men.women.with.opt$se[3,]
	, out.h2.men.women.withLeg.opt$se[3,])
lb_women <- c(out.h2.men.women.non.opt$ci[3,1], out.h2.men.women.with.opt$ci[3,1]
	, out.h2.men.women.withLeg.opt$ci[3,1])
ub_women <- c(out.h2.men.women.non.opt$ci[3,2], out.h2.men.women.with.opt$ci[3,2]
	, out.h2.men.women.withLeg.opt$ci[3,2])

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.h2.men.men.non.opt$coef[3,], beta2 = out.h2.men.women.non.opt$coef[3,]
	, se1 = out.h2.men.men.non.opt$se[3,], se2 = out.h2.men.women.non.opt$se[3,])
zOutStrong <- zDif(beta1 = out.h2.men.men.with.opt$coef[3,], beta2 = out.h2.men.women.with.opt$coef[3,]
	, se1 = out.h2.men.men.with.opt$se[3,], se2 = out.h2.men.women.with.opt$se[3,])
zOutModerate <- zDif(beta1 = out.h2.men.men.withLeg.opt$coef[3,], beta2 = out.h2.men.women.withLeg.opt$coef[3,]
	, se1 = out.h2.men.men.withLeg.opt$se[3,], se2 = out.h2.men.women.withLeg.opt$se[3,])

# Start the plot:
# FIGURE I.63 (a)
pdf('./AppendixIResults/figI63_a.pdf', width = 14)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
	pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-6, 6),
	ylab = "RD Effect of Co-Partisan Mayor",
	xlab = "Electoral Rule Regime", 
	cex.lab = 2)
sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
	y = c(lb_men[i], ub_men[i])))
sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
	y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# Moderate
segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

# Add points for women's coefficients
points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOutWeak$difference, 3), y = betas_men[1] + 3.75, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
, y = betas_men[1] + 2.9, x = 1, cex = 1.5)

text(round(zOutStrong$difference, 3), y = betas_men[2] + 3.95, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
, y = betas_men[2] + 3.0, x = 2, cex = 1.5)

text(round(zOutModerate$difference, 3), y = betas_men[3] + 4.75, x = 3, cex = 1.5)
text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
, y = betas_men[3] + 4, x = 3, cex = 1.5)

# Add legend
legend(x = 0.7, y = -4.5, 
	legend = c('Women', 'Men'), 
	col = c('gray46', 'black'), 
	pch = c(17, 20), 
	lty = c(1, 1),
	bty = 'n',
	cex = 1.5)
dev.off()


################################
######### mayor = women ########
################################
# Weak
out.h2.women.men.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$avg_male_DF_nolar
		, x = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$run_var_fem
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$codeyear
		, all = TRUE)
out.h2.women.women.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$avg_fem_DF_nolar
		, x = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$run_var_fem
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & !(ANO_ELEICAO %in% c(2008, 2012, 2016)))$codeyear
		, all = TRUE)
summary(out.h2.women.men.non.opt)
summary(out.h2.women.women.non.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.women.men.non.opt$coef[3,], beta2 = out.h2.women.women.non.opt$coef[3,]
	, se1 = out.h2.women.men.non.opt$se[3,], se2 = out.h2.women.women.non.opt$se[3,])
outZ
(pnorm(outZ$zscore)) * 2

# Strong
out.h2.women.men.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$avg_male_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$run_var_fem
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$codeyear
		, all = TRUE)
out.h2.women.women.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$avg_fem_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$run_var_fem
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2008, 2012)))$codeyear
		, all = TRUE)
summary(out.h2.women.men.with.opt)
summary(out.h2.women.women.with.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.women.men.with.opt$coef[3,], beta2 = out.h2.women.women.with.opt$coef[3,]
	, se1 = out.h2.women.men.with.opt$se[3,], se2 = out.h2.women.women.with.opt$se[3,])
outZ
(1-pnorm(outZ$zscore)) * 2


# Moderate
out.h2.women.men.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$avg_male_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$run_var_fem
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$codeyear
		, all = TRUE)
out.h2.women.women.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$avg_fem_DF_nolar
		, x = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$run_var_fem
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties & (ANO_ELEICAO %in% c(2016)))$codeyear
		, all = TRUE)
summary(out.h2.women.men.withLeg.opt)
summary(out.h2.women.women.withLeg.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.women.men.withLeg.opt$coef[3,], beta2 = out.h2.women.women.withLeg.opt$coef[3,]
	, se1 = out.h2.women.men.withLeg.opt$se[3,], se2 = out.h2.women.women.withLeg.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2


# Get point estimates, se, and ci for men
betas_men <- c(out.h2.women.men.non.opt$coef[3,], out.h2.women.men.with.opt$coef[3,]
	, out.h2.women.men.withLeg.opt$coef[3,])
se_men <- c(out.h2.women.men.non.opt$se[3,], out.h2.women.men.with.opt$se[3,]
	, out.h2.women.men.withLeg.opt$se[3,])
lb_men <- c(out.h2.women.men.non.opt$ci[3,1], out.h2.women.men.with.opt$ci[3,1]
	, out.h2.women.men.withLeg.opt$ci[3,1])
ub_men <- c(out.h2.women.men.non.opt$ci[3,2], out.h2.women.men.with.opt$ci[3,2]
	, out.h2.women.men.withLeg.opt$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.h2.women.women.non.opt$coef[3,], out.h2.women.women.with.opt$coef[3,]
	, out.h2.women.women.withLeg.opt$coef[3,])
se_women <- c(out.h2.women.women.non.opt$se[3,], out.h2.women.women.with.opt$se[3,]
	, out.h2.women.women.withLeg.opt$se[3,])
lb_women <- c(out.h2.women.women.non.opt$ci[3,1], out.h2.women.women.with.opt$ci[3,1]
	, out.h2.women.women.withLeg.opt$ci[3,1])
ub_women <- c(out.h2.women.women.non.opt$ci[3,2], out.h2.women.women.with.opt$ci[3,2]
	, out.h2.women.women.withLeg.opt$ci[3,2])


# Calculate Differences 
zOutWeak <- zDif(beta1 = out.h2.women.men.non.opt$coef[3,], beta2 = out.h2.women.women.non.opt$coef[3,]
	, se1 = out.h2.women.men.non.opt$se[3,], se2 = out.h2.women.women.non.opt$se[3,])
zOutStrong <- zDif(beta1 = out.h2.women.men.with.opt$coef[3,], beta2 = out.h2.women.women.with.opt$coef[3,]
	, se1 = out.h2.women.men.with.opt$se[3,], se2 = out.h2.women.women.with.opt$se[3,])
zOutModerate <- zDif(beta1 = out.h2.women.men.withLeg.opt$coef[3,], beta2 = out.h2.women.women.withLeg.opt$coef[3,]
	, se1 = out.h2.women.men.withLeg.opt$se[3,], se2 = out.h2.women.women.withLeg.opt$se[3,])

# Start the plot:
# FIGURE I.63 (b)
pdf('./AppendixIResults/figI63_b.pdf', width = 14)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
	pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-6, 6),
	ylab = "RD Effect of Co-Partisan Mayor",
	xlab = "Electoral Rule Regime", 
	cex.lab = 2)
sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
	y = c(lb_men[i], ub_men[i])))
sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
	y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# Moderate
segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

# Add points for women's coefficients
points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOutWeak$difference, 3), y = betas_men[1] + 6, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
, y = betas_men[1] + 5.1, x = 1, cex = 1.5)

text(round(zOutStrong$difference, 3), y = betas_men[2] + 3.95, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
, y = betas_men[2] + 3.05, x = 2, cex = 1.5)

text(round(zOutModerate$difference, 3), y = betas_men[3] + 3.75, x = 3, cex = 1.5)
text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
, y = betas_men[3] + 3, x = 3, cex = 1.5)

# Add legend
legend(x = 0.7, y = -4.5, 
	legend = c('Women', 'Men'), 
	col = c('gray46', 'black'), 
	pch = c(17, 20), 
	lty = c(1, 1),
	bty = 'n',
	cex = 1.5)
dev.off()
